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SUMMARY 


The traditional approach for combustor development, based on cut-and-try 
experimentation, is very costly, and the designer of supersonic combustors 
has encountered difficulties in the acquisition of meaningful data. Recent 
advances in computational capabilities have made it possible to predict complex 
flow fields in combustors. Such computational approaches could supplement the 
traditional experimentation and help overcome some of the difficulties in super- 
sonic combustor development. 

Two parabolic flow computer programs, one based on a finite-difference 
method and the other on a finite-element method, are used at the NASA Langley 
Research Center. Both programs are capable of predicting three-dimensional 
turbulent reacting flow fields in supersonic combustors. A three-dimensional 
turbulent mixing experiment has been conducted for detailed evaluation of the 
two computer programs. To alleviate the difficulties associated with high- 
temperature measurements, a cold (nonreacting) flow experiment was performed to 
study the mixing of helium jets with a supersonic air stream in a rectangular 
duct. Surveys of the flow field were made at two stations downstream of the 
injectors. The flow surveys at the upstream station were used as initial data 
by both computer programs. The surveys at the downstream station provided an 
experimental comparison to assess the relative accuracies of the two programs. 
Computational efficiencies were also compared. 

In this paper, the theoretical foundations of the two computer programs 
are first introduced. The three-dimensional mixing experiment and the applica- 
tion of two computer programs to that experiment are then described. Finally, 
comparisons between the two computations and between the computations and the 
experiment are made. In general, both computer programs predicted the experi- 
mental results and data trends reasonably well. However, the comparison indi- 
cated that the finite-difference program was more accurate in computation and 
more efficient in both computer storage and computing time than the present 
version of the finite-element program. 


INTRODUCTION 

Recent advances in computational techniques and computer capabilities have 
made it possible to predict three-dimensional turbulent reacting flow fields. 
Such analytical predictions are extremely valuable in the development of super- 
sonic combustors for planned hypersonic airbreathing engines. The traditional 
approach for such combustor development has been based on expensive cut-and-try 
exper imentation . 

At present, two computer programs, SHIP (ref. 1) and CCMOC (ref. 2), which 
are capable of predicting three-dimensional turbulent flow fields in supersonic 
combustors, are operational at the NASA Langley Research Center. SHIP is based 
on a finite-difference algorithm of Spalding (ref. 3) and COMOC is based on a 


finite-element algorithm of Baker (ref. 4). The mathematical foundations and 
the relative merits of the two algorithms are well known; successful computa- 
tions using these algorithms have appeared in the open literature (for example, 
refs. 2, 3, and 5 to 8) . However, direct comparisons between these two com- 
puter programs have not been made to assess their computational capabilities, 
turbulence models, numerical schemes, and results. The objective of this paper 
is to make such direct comparisons to assess both programs as combustor design 
tools. To aid in accomplishing this objective, a three-dimensional mixing flow 
experiment has been conducted and both computer programs have been applied to 
predict the experimental flow field. 

In the past, the two computer programs were applied separately to a mix- 
ing flow field with normal injection of hydrogen into a supersonic airstream 
(ref. 9). However, meaningful comparisons were not possible because of differ- 
ent approaches to handling the same flow problem. The finite-element program 
was applied downstream of the injection by modeling the normal injection with 
an equivalent one-dimensional virtual source (ref. 7) . Since this computation 
depends greatly on the appropriateness of one-dimensional modeling and the 
accuracy of such modeling is difficult to estimate, comparison of such a compu- 
tation with experimental data can hardly establish the capability of the com- 
puter program. On the other hand, the finite-difference program was applied 
across the normal injection region even though there was a limitation due to 
the parabolic flow assumption (ref. 10) . It was later found that such a compu- 
tation could yield inconsistent results because of the presence of flow recir- 
culation. By using the near-field experimental measurements of reference 9 as 
input data, the same program was recently applied to the flow field downstream 
of the injection and the recirculating regions (ref. 11) . Because of the insuf- 
ficiently detailed initial data, only qualitative agreement between the compu- 
tation and experiment was obtained. Therefore, to evaluate and establish the 
three-dimensional capabilities of the present computer programs, well-designed 
experiments are required which provide detailed flow-field surveys at several 
stations in a supersonic combustor. 

It is well recognized that there are many technical difficulties associated 
with experimental measurements in the supersonic turbulent reacting flow fields 
of a combustor. Because of the highly turbulent, high-temperature environment, 
the acquisition of meaningful data is extremely difficult. To alleviate the 
difficulties associated with high temperature, a three-dimensional cold (nonre- 
acting) flow experiment was performed to study the mixing of helium jets with 
the supersonic airstream in a rectangular duct. Surveys of the flow field were 
made at two stations downstream of the injectors, one close to the injectors 
and the other further downstream. The flow surveys at the upstream station 
were then used as initial data in both computer programs. The surveys at the 
downstream station provided an experimental comparison to assess the relative 
accuracies and computational efficiencies of the computer programs. 

The mathematical foundations and numerical schemes of the two computer 
programs are presented in the first section. Details from earlier work 
(refs. 1 to 4) have been collected and expanded to provide a convenient com- 
parison of the theories. The three-dimensional mixing experiment is then 
described relative to test conditions, procedure, data acquisition, and data 
reduction. Applications of the two computer programs to the experiment are 
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II 


presented, and comparisons of flow-field quantities are made between the two 
computations and between the computations and the experiment. Finally, the 
capabilities of the two programs to compute three-dimensional mixing flow 
fields are discussed and conclusions are drawn. 

SYMBOLS 

A duct cross-sectional area 

Ajj computational cross-sectional area 

A}], Ag, A£,Af^ V coefficients in difference equation ( 5 ) 

® J 

a ( 1 ) ,a ^ 2 ) ,a boundary-condition coefficients 
ad)' = ad)/a(2) 

a(3)' = a(3)/a(2) 

a,b,c functional coefficients 

0^,02, Ci3,C]^ empirical constants 

D Van Driest damping factor 

f mass fraction; function of dependent or independent variables 

g function of dependent or independent variables 

H total enthalpy 

h static enthalpy 

K general diffusion coefficient 

k turbulence kinetic energy 

L,)l differential operators 

2 jg mixing length 

M Mach number 

m mass flow rate 

Np^ Prandtl number 

ny,n2 unit normal to y- and z-direction, respectively 
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p 

p 

Q 

q 

R 

3r 

S 

T 


u,v,w 


V 
W 

XfY»z 

r 

y 

«b 

e 

% 

ic 

X 

y 

V 

p 

T 


node point 
pressure 

nodal value vector of dependent variable 

general dependent flow variable 

solution domain 

solution domain boundary 

mth finite-element subdomain 

mth subdomain boundary 

source term; surface area 

temperature 

velocity component in x-, y-, and z-direction, respectively 
volume 

weighting function 

rectangular coordinates 

exchange coefficient 

ratio of specific heats 

boundary-layer thickness 

mixing-layer thickness 

turbulence dissipation energy rate 

mixing efficiency 

coefficient in equation (7a) 

algebraic multiplier; empirical constant 

viscosity 

kinematic viscosity 

density 

shear stress 
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general functions 

4> 

space functional 

Subscr ipts : 

A 

air 

av 

average 

eff 

effective 

He 

helium 

I 

injectant 

IL 

laminar 

m 

finite-element index 

n 

node index 

T 

total 

t 

turbulent 

w 

wall 

Superscripts: 

a 

species a 

* 

approximate solution 


restrained to boundary of computational domain 


DESCRIPTION OF NUMERICAL METHODS 

In this section the theoretical bases and numerical schemes of the finite- 
difference (SHIP) and finite-element (COMOC) computer programs are described. 
Both programs are developed on the basis of an Euler ian formulation in a rec- 
tangular coordinate system (x, y, and z) with the x-axis in the main flow 
direction. The mean flow velocity components (u, v, and w) , pressure (p) , 
total enthalpy (H ) , and mass fraction (f ) of a three-dimensional turbulent 
mixing flow field are governed approximately by the Navier-Stokes equations 
together with a species equation. Both programs can consider flow with or 
without chemical reactions. Moreover, the flow is assumed to be composed of 
perfect gases with specific heats that are functions of temperature and species. 


The effect of turbulence is introduced by replacing the laminar vis- 
cosity (Ug,) or the laminar exchange coefficient by an effective viscosity 
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(V^ff = iJj. + yjjj or the corresponding effective exchange coefficient. The tur- 
bulent viscosity (P^) is determined by means of a turbulence model. In the two 
computer programs, different turbulence models are used at the present time and 
they will be discussed separately. 

A parabolic flow assumption is used in both programs to simplify their 
formulations. Physically, this assumption is valid for a flow when there 
exists a predominant flow direction, when the diffusions of mass, momentum, 
energy, etc. , in this direction are negligible compared with the corresponding 
convections, and vAien the downstream pressure has little effect on the upstream 
flow field. Mathematically, the set of governing equations reduces to a system 
of parabolic-type equations. Numerically, these equations can be solved in 
succeeding cross-stream (y-z) planes progressing in the main flow direction. 
Thus, a three-dimensional problem requires only two-dimensional computer stor- 
age, and computer time and storage are greatly reduced. However, because of 
the parabolic flow assumption, the range of application of the programs is also 
limited accordingly. 

The different features of two programs are described separately in the 
following discussion. 


Finite-Difference Computer Program (SHIP) 

The flow field considered in the finite-difference computer program (SHIP) 
is defined by a rectangular parallelepiped. Any of the four lateral boundaries 
can be a free, symmetry, or wall surface; for walls, the boundaries are allowed 
to expand or contract smoothly in the main flow direction. The main flow can 
be either subsonic or supersonic. 

The governing equations for three-dimensional parabolic flow can be writ- 
ten as follows (see also ref. 11) : 

Continuity 


9 3 8 

■^(Pu) + -^(Pv) + — (Pw) 

dx oy dz 


= 0 


(la) 


x-momentum 


a ^ a a 

— (Pu2) + — (pvu) + — (Pwu) 

ax ay az 


ax 




(lb) 


y-momentum 

a a a 

— (Puv) + — (Pv2) + — (Pwv) 

ox oy dz 


ap a 

2 / av 

aw\ 

a 

/av aw\ 

ay ay 

I ^ - 

w 

+ — 
az 

^ Ty) 

•m — 


(Ic) 
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z-momentura 


3 3 3 

— (puw) + — (pvw) + — (PW‘‘) 

3x 3y 3z 


3p 

3 

/ 

'3v 

3w 

+ 

V»eff( 

+ 

3z 

3y 

^3z 

3y 


+ 


3 

3z 


2 

~ l^eff 


3w 

3z 



(Id) 


Energy 


3 3 3 

— (puH) + — (pvH) + — (pwH) 

3x 3y 3z 


Meff 


3h^ 


V*ef f 


3H) 


9y\Npr,eff,H \Npi-, eff , H 


+ 


3 

3y| 


V^eff 


N 


Pr,eff ,H 


■(Npr,eff,H 




+ 


3 

3z 


^^eff 


Npr,eff,H 


(Npr,eff,H 


3 /u^ + v2 

D— 

3z\ 2 



+ Vieff 



+ 


/3w\^ 3v 3w 

\3z/ 3z 3y 


(le) 


Species^ 


3 3 3 

— (puf) + — (pvf) + — (puf) 
3x 3y 3z 


^ef f 


3f 


9y\Npr,eff,f 3yy 


V*ef f 


3z\Npr,eff,f 



(If) 


Here, **Pr,eff,H Npr,eff,f ®re the effective Prandtl numbers for total 

enthalpy and mass fraction, respectively. The general effective exchange coef- 
fient r is expressed in the following form. 


Peff 

T = (2) 

Npr,eff **pr,Jl %r,t 

with Npr^i!, and Npj-^t being the laminar and turbulent Prandtl numbers, respec- 
tively. Since Npr,£ and Npr,t of bbe order of unity and 

fully developed turbulent flow, T ~ yt/^Pr,t* Note that some terms in the set 
of equations (1) which were omitted in reference 1 have been included for com- 
pleteness, and the finite-difference program has also been updated accordingly. 


^Only one differential equation for the mass fraction is solved; remaining 
species are calculated on the basis of stoichiometric considerations. 
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A "k-e" two-equation turbulence model is used (ref. 1) . From a dimensional 
analysis, the turbulent viscosity is proportional to the flow density, a 

turbulence velocity scale (e.g., the square root of the turbulence kinetic 
energy k) , and a turbulence length scale 1 . In high Reynolds number flows, 

I • k^/2/e with e being the turbulence dissipation energy rate; hence, 

Ut = CijPk^/e where Cjj is an empirical constant and k and e are governed 
by a set of transport equations. Under the parabolic flow assumption, these 
transport equations are 


9 9 9 9 / ^eff 3k\ 3 / *^eff 3k' 

— (Puk) + — (Pvk) + — (Pwk) = — + — 

3x 9y 3z 9y\Npr,eff,k ^ 3z\Npr,eff,k 


+ Pt\ 2 


'9v\ /3w' 

,3y/ \3zy 


/3u\ / 3u 

+ — + — 

\9y/ \9z, 


/ 3v 3w' 


(3a) 


9 9 9 9 / ^‘eff 3e\ 3 / ^eff 3e' 

— (Pue) + — (Pv£) + — (Pwe) = — — + — 

3x 3y 9z 3y\Npr,eff,e ^ 5z\Npr,eff,e 


+ Cl - PtS 2 
k 


9v ' 

9y/ 


3w' 
3z y 




\ h) 


/3u\^ / 3v 3w\^ I pe2 


(3b) 


where Ci and C 2 are empirical constants associated with the k-e two- 
equation turbulence model. 

The set of simplified parabolic flow equations can be cast in a conserva- 
tive form: 


9 

— (Puq) 
3x 



S 


q 


8 


(4) 



%^ere q is a general dependent variable and and Sq are the exchange 

coefficient and source term, respectively, associated with the dependent vari- 
able q. nhen q equals 1, u, v, w, H, f, k, or e, equation (4) cor- 

responds to the equations for continuity, the three components of momentum, 
energy, species, turbulence kinetic energy, or turbulence dissipation energy 
rate. Numerical computation is based on a finite-difference formulation of 
equation (4). A "staggered" grid system is used in the cross-stream (y-z) 
plane and variable grid spacings are allowed. Figure 1 shows such a staggered 
grid system. In the figure, P indicates an arbitrary node point; its four 
neighboring points are denoted by E, W, N, and S for east, west, north, 

and south. All variables at an arbitrary node are stored at the location P 

except the transverse velocity components v and w: v is stored at the mid- 

point between W and P, and w is stored at the midpoint between S and P, 
as shown by the arrows in figure 1. Thus, control volumes at each node point 
are different for v euid w compared with those for the other dependent vari- 
ables. By taking volume integrations of equation (4) over respective control 
volumes, a set of difference equations can be obtained. 

3 

The volume integrations of the terms — (Puq) and So in equation (4) 

3x ^ 

are performed by assuming that the values of q, p, and u at a node point P 
are constant over the entire control volume. The term 3/3x contains the 
difference of values of stations x and x + Ax. The volume integrations of 
the other two terms in equation (4) give rise to the surface integrals of the 
convective and diffusive fluxes across the boundaries of the control volume. 
Proper representation of these terms is essential to the convergence of the 
numerical computation. To provide numerical convergence and accuracy, a 
"hybrid" scheme is used (ref. 12) which is a combination of central and 
upwind differences, when the integrations of various terms in equation (4) 
are expressed in the manner just described, the general form of a difference 
equation at an arbitrary node point P is written in the following form: 

*3p = ^‘IN Ajiqg + Aff<^ + B (5) 

where the values of the q's pertain to station x + Ax and the values of the 
A's and B pertain to station x. The subscripts N, S, E, and W denote 
the neighboring north, south, east, and west nodes (fig. 1). 

The set of difference equations (eq. (5)) together with other auxiliary 
relations are solved by a so-called SIMPLE (£emi-mplicit method for £ressure 
jinked equations) procedure (ref. 1) . Although the governing equations are 
highly nonlinear, an economical noniterative marching procedure in the stream- 
wise direction is followed. The three velocity components are solved first 
from their respective equations (i.e., q = u, v, and w) in terms of a 
guessed pressure field. Then a pressure correction is obtained from an equa- 
tion derived from the continuity equation and having a form similar to equa- 
tion (5) . After the pressure field eind the three velocity components have been 
corrected, the difference equations for H, f, k, and e are solved sequen- 
tially. Temperature, density, and mass fraction for each species and other 
auxiliary variables are determined noniteratively. 
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In the finite-difference program, the boundary conditions are specified by 
the values of appropriate fluxes across the boundaries. For a free-strecim or 
symmetry boundary, fluxes are set automatically to zero. To avoid using a large 
number of grid points for calculating large gradients near a solid wall, the 
boundary conditions are described by the law of the wall at near-boundary grid 
points (ref. 13). Once the boundary conditions are properly specified, each of 
the difference equations (eq. (5)) is solved in a line-by-line iterative fash- 
ion in the cross-stream (y-z) plane by the successive application of a standard 
tridiagonal matrix algorithm in the y- and z-directions. 

The accuracy of the computation is determined by examining the conserva- 
tion of mass flow in each cross-stream plane and the mass residue at each node 
point; both the mass residue at each node point and the mass-flow Imbalance in 
the entire cross-stream plane are expected to be small. The accuracy may be 
improved by increasing the number of grid points in the cross-stream plane, 
increasing the number of iterations for solving difference equations, or using 
smaller forward steps in the streamwise direction. 


Finite-Element Computer Program (COMOC) 

The finite-element computer program (COMCX:) considers flow in a constant- 
area rectangular duct with arbitrary boundary conditions on either the four 
duct walls, symmetry planes lying within the duct, or any combination of these 
boundaries. The flow can be either subsonic or supersonic, with the only 
restriction being that streamwise diffusion is negligible so that a parabolic 
character is maintained. 

Either a mixing-length or a k-e two-equation turbulence model (ref. 2) 
is available in COM(X: to calculate the effective turbulent viscosity. The k-e 
model requires solving two partial differential equations in addition to the 
governing equations. Even though the increase in the number of equations is 
small, the equations for the turbulence kinetic energy and dissipation energy 
rate are quite "stiff" for the COMOC integration algorithm. Although an effort 
is currently underway to streamline the integration technique, the present work 
used the turbulence model based on mixing-length theory (ref. 2) to improve 
computational efficiency for the large domains being considered. 

The parabolic Navier-Stokes, energy, and species equations describing the 
flow are given in reference 2. In the notation of the present paper, they are 

Continuity 


9 3 3 

— (Pu) + — (Pv) + — (Pw) = 0 (6a) 

ox 3y 3z 


x-momentum 


P 


3u 


3u 

+ V — 

3y 


+ 






(6b) 
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y-momentum 


I 3v 9v 9v\ 

p u — + V — + w — = 

y 9x 9y 9z/ 


9p 9 / 9v'' 


9 / 9v' 

9y ^ 9yl^^^^ 9yJ ^ 9zy 


(6c) 


z-momentum 


/ 9w 9w 9w ' 

p u — + V — + w — 
y 9x 9y 9z/ 


JP . 3 / 

9z 




9w' 


^ 3yy ^ ^y) * 9z\^^®^^ 9z^ 


(6d) 


Energy 


/ 9h 9H 9h 

pu — + v — + w — 

y 9x 9y 9zy 


9 / V^eff 9 h\ 9 / ^eff 9 h' 

3y\Npr,eff,H 3y/ SzlNpr^eff ,H 3z 


9y 


1 ~ )%ff 3 o o o 

(uZ + 

Npr,eff,H 2 9y 


3 

9z 


1 - Npr,eff,H Meff 3 . . . 

(u^ + v2 + w2) 

%r,eff,H 2 9z 


9 fNpj-,eff,f - Npr,eff,H ^ 9f “ ' 

— Peff 21 

3yy Npr,eff,fNpr,eff,H a 3y 


9 /Npr,eff,f - Npr,eff,H 
3zy Npr,eff,fNpr,eff,H 


Peff 


h« 


ot 


9f“ 

9z 


(6e) 


Species 
9f» 


p u 


9x 


9ftt 3fa> 

+ V + w 


Ueff 9fa\ 9 
+ 


Peff 9f“' 


9y 9z j 3y\Npr,eff,f 3y / 3z\Npr,eff,f 3z 


+ S“ 


(6f) 


where f®, h“, and S*^ are the mass fraction, static enthalpy, and source term, 

respectively, of species a (for nonreacting flows, 8°^ = 0) . The effective 
Prandtl number for the mass fraction Npj-^eff^f (the effective Schmidt number) 
is assumed to be the same for all species. Equations (6b) to (6e) do not form a 
completely parabolic system because certain cross-derivative shear terms in the 
y- and z-momentum equations and viscous dissipation terms in the energy equation 
have not been considered. These terms were assumed to be small in the develop- 
ment of the finite-element program and therefore were dropped. They are cur- 
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rently being added/ however, in a newer version of the program. Equation (6) 
can be expressed in general form as 


L(q) 



/ 3q 3q 3q \ 


g(q,x) = 0 


(7a) 


with generalized boundary conditions 


^(q) 


a(l)q(x,y,z) 


+ 


a(2)K(q) 



q(x,y,z)ny 


3 

+ — q(x,y,z)n2 

dz 


a (3) =0 (7b) 


where q is any dependent variable, K is a generalized diffusion coefficient, 
f and g are functions of the specified variables, and <, a(^), a(2), and 

a (3) are specifiable constants. The superscript bar restrains y and z to 
the boundaries of the computational domain. 


A computational domain is defined identically with the duct cross-sectional 
area at the initial station of the solution domain. This area is then discre- 
tized, as shown in figure 2, with triangular finite elements sized by the user 
with respect to the initial rate of change of the variables q in each element 
(solution subdomain) . The triangles are chosen small for high resolution where 
the gradients of q are large; they are chosen large for computational economy 
where the gradients of q are small. This finite-element scheme provides a 
convenient means of transforming the partial differential equations describing 
the system into coupled ordinary differential equations that can be more readily 
solved. Within each triangular element, the variables q are assumed to vary 
in a linear fashion described by 

qm(x,y,z) = am + bmy + CmZ (8) 

where qm is the approximate solution to q in the mth element and a, b, 
and c define the values of q^ at the three corners of the triangle. For 
subsequent use, equation (8) is more conveniently expressed in vector form as 

qm(x,y,z) = {4> (y, z) }'^{Qm (x) > (9a) 

where 

,z) } = 




{Qm(x)> = 



(9c) 


euid superscript T indicates the transpose of the vector. 
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To determine the Q's at each node, which are the solutions sought, the 
approximate solutions q* are substituted into the differential equations 
(eq. (7a)) constrained by boundary conditions (eq. (7b)). If the values of 
L(q*) and ^>(q*), referred to as residuals, can be reduced collectively to 
zero over the solution subdomain, the approximate solution then approaches the 
exact solution of equation (7a) , as the element size is reduced to zero in the 
limit. 

Values of the residuals are best minimized 
element subdomain with weight W and requiring 
is, 

W L(qm) dV - X (P W il (q^) dS 

Rm 

This procedure is commonly referred to as the method of weighted residuals. 
Here, denotes the element subdomain, the subdomain boundary, 3R 

the computational domain boundary, and n the intersection of these two bound- 
aries. The second integration is performed only when an element lies on the 
boundaries of the duct cross section being computed. Weighting allows the 
functional dependence of q* to be felt throughout the integral domain. The 
unknown algebraic multiplier X has been included (ref. 2) for use later in 
simplifying the equation system. 

In order to proceed with equation (10), appropriate weights on L and Z 
must be defined. The Galerkin procedure, applied here, chooses the weights to 
be identical with (j) (eq. (9b)). Experience has shown that this procedure 
yields the most appropriate weighting to equation (10) consistent with the sys- 
tem being solved. Equation (10) then becomes 

{(|)(y,z)} L(q*) dV - \ (b {4)(y,z)} S, (q*) dS = {0} (11) 

'^Rm *^9Rn,n3R 

Equation (11) describes the variation of the dependent variables q over 
each finite element within the computational domain. To describe q over 
the entire domain, equation (11) with appropriate simplifying modifications 
(refs. 2, 7, and 14) must be summed (symbol U) over all the finite elements. 
The details of this development are rather involved and are presented in the 
appendix; the result is 


by integrating over the entire 
that the result vanish; that 

= 0 (10) 




3{(j)} 9qm 

K + 

. 37 3y 


‘^3R„in3R 


3{(|)} 

K 

3z 


3z 


dV + 





( 12 ) 
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Equation (12) is the basis of the solution algorithm for all dependent variables 
of the COVOC finite-element computer program. It is applied directly to the 
three momentum equations, the energy equation, and the species equation. A sys- 
tem of simultaneous first-order ordinary differential equations in the x 
(streamwise) coordinate direction then results for the dependent variables. 

These equations are solved at each streamwise station using a predictor- 
corrector numerical integration scheme (ref. 7). The value of each dependent 
variable is predicted at the first downstream station (first step) using a sim- 
ple Euler scheme based on the variable's streamwise derivative at the initial 
station. Predictions at subsequent steps are obtained from the predictor for- 
mula. At each step, the predictions are then corrected by the more accurate 
corrector formula. The scheme is absolutely stable and has an error of the 
order of the square of the step size in the streamwise direction. The step 
size is adjusted with respect to a maximum relative truncation error and is 
continually increased up to this limit while downstream stepping to reduce the 
required computer time. When the error limit is exceeded, the preceding step 
size to a new station is reduced and the values of the dependent variables at 
the new station are recalculated. 

As the solution marches downstream, continuity is enforced by insuring 
that mass is conserved at each streamwise station. Conservation is guaranteed 
by correcting the streamwise pressure gradient at each station, so that the 
computational cross-sectional area of the flow A^, matches as closely as pos- 
sible the actual cross-sectional area A of the duct. This computational 
cross-sectional area is directly related to the total mass flow rate m by 


Ac = 



(13) 


where Pj, and Uj^ are local values of density and streamwise velocity at each 
of the N nodes forming the finite-element grid. The change in computational 
area necessary to match the actual area is then. 


dA^ A A^j 
dx Ax 


(14) 


where Ax is the step size. The derivative dAj,/dx is directly related to 
the streamwise pressure gradient dp/dx by (ref. 13) 
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where is the average wall shear over the step, Tgy is the area-averaged 

temperature, is the area-averaged Mach number over the ‘cross section, and 

Y is the ratio of specific heats. COMOC currently considers only this stream- 
wise variation of pressure. However, the capability to calculate both trans- 
verse pressure fields is now being implemented in the program. 

Specifiable boundary conditions in COMCX: are of three basic types. The 
symmetry (vanishing gradient) boundary condition, where the normal derivative 
of any dependent variable vanishes, is the default boundary condition in COMOC 
(see appendix) . Alternately, the value of dependent variables can be held at 
its initial value along any boundary. Finally, a finite value of the normal 
derivative of a dependent variable along any boundary can be specified. This 
"derivative boundary condition" is useful in turbulent flows and will be 
described more fully later. 


EXPERIMENTAL TEST CASE 
Apparatus and Test Conditions 

A three-dimensional turbulent mixing experiment has been conducted to pro- 
vide detailed data for comparison of the two computer programs. The experiment 
involves the mixing of helium jets with a supersonic airstream in a rectangular 
combustor duct. A sketch of the injector strut and combustor duct is presented 
in figure 3. The combustor duct has a constant cross section 0.0381 m by 
0.17 m. The first section of the duct contains the injection strut, which is 
centered in the shorter dimension and spans the longer dimension of the duct 
cross section. The combustor duct incorporates hinge joints at the 0.309-m 
station for making adjustments to the combustor geometry. Flow to the combus- 
tor duct is provided by a rectangular Mach 2.99 nozzle at a total temperature 
of 294 K and a total pressure of 1.6 MPa (see table I). 

Details of the injector strut are presented in figure 4. Previous tests 
of this strut, reported in reference 15, were performed with an area ratio of 2 
over the aft 0.479-m section of the combustor. This strut was designed for 
combustion tests, and therefore, the leading edge has provisions for water 
cooling. The strut consists of a 6° half -angle wedge with a maximum thickness 
of 0.01 m and a 0. 0016-m-radius leading edge. The aft body has constant thick- 
ness and a blunt base. Five equally spaced fuel injectors located on the base 
divide the flow into five nearly square (0.0340 m by 0.0381 m) regions. Four 
of the injectors are 10° half -angle conical nozzles which operate slightly 
overexpanded. The other injector has two sonic jets directed toward each other 
with a combined jet cross-sectional area equal to the throat area of the other 
jets. The injectant is supplied through the passages labeled "Helium" in 
figure 4. As listed in table I, helium is injected at a total temperature of 
294 K and a total pressure of 3.27 or 3.79 MPa, which produce a bulk helium-to- 
air mass-flow ratio of 0.024 or 0.027, respectively. 
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Instrumentation 


Measurements included total injectant flow rate, total pressure and tem- 
perature of air and injectant, combustor wall pressures, and instream surveys 
of static pressure, pitot pressure, and gas composition. Comparisons between 
the injectant flow rate measured by a venturi in the helium supply line and the 
flow rate calculated by one-dimensional isentropic compressible flow relations 
at the injector throat show that the bulk discharge coefficient for the five 
jets was approximately 0.81. Combustor wall pressures were measured along the 
entire length of the upper wall, primarily on the wall center line. Instream 
flow surveys of pitot pressure and gas composition were made with a 13-probe 
rake which incorporated internal expansion pitot probes, as illustrated in fig- 
ure 5(a). The probes were mounted at a center-to-center spacing of 0.00635 m. 
The pitot probe utilizes an internal expansion design to aid in collecting 
accurate gas samples by reducing or eliminating tip spillage. All 13 pitot 
pressures were recorded simultaneously; then 13 gas sample bottles were filled. 
These bottles were analyzed on a process gas chromatograph between survey runs. 
Instream static pressure was obtained with a probe rake incorporating the static 
probe tips illustrated in figure 5(b). This static probe tip is less suscepti- 
ble to errors from misalignment with the flow direction than conventional static 
probes (ref. 16) . Static tips were located on the rake at a center-to-center 
spacing of 0.0127 m (double the pitot spacing) to eliminate flow interference 
between the tips. Thus, two passes were required to complete static pressure 
flow mapping at each survey station. 

The metered air and helium mass flow rates are estimated to have an accu- 
racy within ±2 percent. Accuracies of pitot pressure, static pressure, and gas 
concentration are estimated to be within ±2 percent, ±4 percent, and ±1 percent, 
respectively. Another source of error in the data acquisition is repeatability 
of the test conditions for different surveys. The air and helium mass flows 
varied by about ±5 percent from test to test, so the measured pitot and static 
pressures were adjusted by the ratio of the actual air total pressure to mean 
air total pressure. 


Experimental Results 

Flow-field surveys were made at two streamwise stations, x = 0.0421 m 
and X = 0.7369 m downstream from the base of the injection strut. At the 
upstream station (x = 0.0421 m) , the surveys were concentrated in the region of 
the flow occupied by the injectant. At the downstream station (x = 0.7369 m) , 
the entire flow field was surveyed. Measured pitot and static pressures and 
gas concentrations were reduced to generate tabulated data, flow-field profiles, 
and contour plots of various flow properties. 

Helium mass fraction contours for the entire duct cross section at the 
downstream station (x = 0.7369 m) are presented in figure 6. Corresponding to 
the five injectors are five almost equally spaced mixing regions with distinct 
boundaries. Dashed vertical lines, located midway between the projected loca- 
tions of the jet center lines, emphasize the mixing regions. Integrations of 
air and helium mass flows in each of the five mixing regions show that the mass 
flows are almost the same in each region, although at the downstream station the 
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jets have already merged with their neighboring jets. As shown in figure 6, 
the flow field downstream of the crossed injector (second from left) is more 
uniform than those of the other parallel injectors. All four parallel-injection 
mixing regions are quite similar with the exception of slight shifts from the 
projected injection location (solid symbols in fig. 6) . The center jet region 
was surveyed more extensively than the other jet regions. Therefore, the cen- 
ter jet was chosen as a typical flow field to be computed by the two computer 
programs. 

The reduced survey data in the center jet region for the streamwise veloc- 
ity, static pressure, temperature, and helium mass fraction are tabulated in 
tables II and III for the upstream (x = 0.0421 m) and the downstream 
(x = 0.7369 m) stations, respectively. These data have been adjusted to the 
bulk flow pressures listed in table I and were reduced by assuming constant 
total temperature. In the following discussions, detailed surveys of the cen- 
ter mixing region are presented. 

Figure 7 Illustrates the center jet mixing region relative to the upstream 
projection of the strut and center fuel injector. The top and bottom boundaries 
are formed by the combustor wall, and the side boundaries by the center planes 
between adjacent jets. Vertical locations of the probe tips and horizontal 
locations of the rake at the upstream survey station (x = 0.0421 m) are repre- 
sented by the tick marks on the respective boundaries. Thus, the intersections 
of horizontal and vertical lines through the tick marks represent the locations 
of survey points. 

Figures 8(a) and 8(b) present the helium mass fraction contours at the 
upstream and downstream survey stations, respectively. At the upstream survey 
station, the five individual jets are still separate, each having spread over 
only about half of the mixing-region width. At the downstream station, the 
flow is nearly stratified in the vertical direction because of merging between 
adjacent jets. The maximum helium mass fraction measured at the upstream sta- 
tion is 0.36; at the downstream station, it is 0.035. 

Helium mass flow rate contours are very similar in appearance to the mass 
fraction contours of figure 8 and are not presented. Helium mass flow rate 
contours were numerically integrated and the result was nondimensionalized by 
the metered mass flow rate from one jet. The result, presented in figure 8 by 
the symbol mjjg, shows a 16- and 10-percent deficiency in helium mass-flow 
measurement at the upstream and downstream stations, respectively. This is 
consistent with previous cold mixing studies where surveys taken close to the 
jet in regions of sharp concentration and velocity gradients exhibit significant 
injectant mass-flow deficiencies. A 12-percent deficiency in air mass-flow 
measurement was observed by integrations of the air mass flow rate contours. 

Velocity contours are presented in figures 9 (a) and 9 (b) at the upstream 
and downstream survey stations, respectively. The helium jet exit velocity is 
1539 m/sec, and the air velocity ahead of the strut is 617 m/sec. At the up- 
stream station surveyed, the jet has a much higher velocity than the surround- 
ing air. The peak velocity in the jet flow is about 1100 m/sec; outside the 
central jet region, the velocity is in the range of 366 to 488 m/sec. This air 
velocity deficiency from 617 m/sec ahead of the strut is indicative of the 
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pressure loss caused by strut- induced shock waves. At the downstream station, 
the velocity is more uniform emd has decayed as a result of viscous losses in 
the duct. The velocity ranges from 335 m/sec near the wall to 457 m/sec at the 
combustor center line euid appears similar to turbulent-boundary-layer-type flow. 

Static pressure contours are presented in figure 10 for the two survey 
stations. At the upstream station (fig. 10(a)), pressure in the region of 
the injectant is lower than in the surrounding flow. The pressure increases 
slightly toward the wall, as expected, because of the wave structure produced 
by the strut. At the downstream station (fig. 10(b)), the static pressure is 
more random but at a higher value because of viscous losses in the duct. 

Static temperature contours, reduced from the survey data by assuming con- 
stant total temperature, are presented in figure 11 for the two survey stations. 
Because the total temperature is assumed to be constant, these surveys are sim- 
ilar to the velocity contours. At the upstream station (fig. 11(a)), the mini- 
mum static temperature at the center of the helium jet is about 139 K. The 
maximum static temperature adjacent to the wall is about 256 K. At the down- 
stream station (fig. 11(b)), the minimum and maximum values are about 194 K and 
244 K. 


NUMERICAL COMPUTATIONS 

In this section, the finite-difference and finite-element computer programs 
are applied to the helium mixing experiment described in the previous section. 
Because of the symmetric nature of the flow field downstream of each jet, only 
the flow field downstream of the central jet is considered in the present calcu- 
lation. In order to take full advantage of the maximum storage in the computer 
programs, the finite-difference program computed only half of the central jet 
(region ABCDEFA in fig. 12), aind the finite-element program computed one-quarter 
of the central jet (region ABCFA in fig. 12). The approximate symmetry of the 
central jet flow field allowed this reduction in the computational domains. 

Since for comparison of the two programs the two computations must be based on 
the same set of initial conditions, the acquisition of the initial conditions 
will be described first. Then, the application of two programs to the problem 
with respect to their grid arrangements, initial and boundary conditions, ini- 
tial estimates of turbulence kinetic energy and dissipation rate, and computer 
time and storage required is described. Finally, computational results of both 
computer programs are presented for comparison with each other and with the 
experimental data. The essential differences between the two computations from 
SHIP and COMOC are presented in table IV. 


Initial Data 

As described in the section, "Experimental Test Case," the measurements 
of static and pitot pressures and helium mass fraction were made at two stream- 
wise stations. From these measurements, the streamwise velocity, temperature, 
pressure, air and helium mass fractions, and density of the mixture were 
deduced. In the center jet region, however, there are only 9x7 points of 
measurement, too coarse to represent the details of the initial conditions for 


18 



the computations. Therefore/ the experimental data (streamwise velocity/ tem- 
perature/ pressure/ and helium mass fraction) were interpolated according to 
the desired grid arrangement in the initial cross-stream plane by using cubic- 
spline interpolation. Initial profiles of streamwise velocity/ temperature/ 
pressure/ cuid helium mass fraction from the center plane (z = 0) to the duct 
wall (z = 0.01905 m) at several stations along the y-axis are shown in fig- 
ure 13. The symbols are the experimental data and the lines are the profiles 
Interpolated from these data. Note that at the initial station/ there were 
shocks; however/ because of the limited measured data/ the positions and inten- 
sities of these shocks cannot be easily determined. Therefore/ the initial 
profiles of the interpolated flow variables show smooth variations over the 
entire cross-stream plane. It is believed that in the parabolic flow field/ 
the error due to such approximations does not produce significant downstream 
effects since any discontinuities will be diffused (smoothed out) within a 
short distance downstream. 


Finite-Difference Program 

In the finite-difference computer program (SHIP ) , the initial conditions 
required are the three velocity components/ temperature/ pressure/ and mass 
fraction of all species. From these flow variables/ the values of density/ 
total enthalpy/ and turbulence kinetic energy and dissipation rate can be 
either calculated or estimated. In the present computation, the two lateral 
velocity components were taken to be initially zero, because they were not mea- 
sured and were assumed to be small compared with the streamwise velocity com- 
ponent. The helium mass flow and the total mass flow computed from the initial 
profiles are 0.01193 kg/sec and 0.5975 kg/sec, respectively, whereas they were 
0.01425 kg/sec and 0.5960 kg/sec in the experiment. 

The initial turbulence kinetic energy and dissipation rate also were not 
measured; hence, they must be estimated. One method, suggested in reference 17, 
is to use an estimated shear-stress profile. For a two-dimensional subsonic 
turbulent boundary layer or turbulent mixing layer of homogeneous medium, it 
was found (ref. 11) that an estimate based on the Prandtl mixing-length hypoth- 
esis is adequate. For a turbulent mixing layer of different mediums, however, 
the mixing length seems to depend on the variation of density (ref. 11) . Since 
there are no rigorous conventional methods available to determine the turbulence 
kinetic energy and dissipation rate from the local flow variables in a three- 
dimensional flow, the following estimate was used in the present finite- 
difference computation. 

It was assumed that the turbulence kinetic energy at the initial sta- 
tion was composed of two parts, one due to the background value and the other 
due to the shear stress, as suggested in reference 17 for two-dimensional 
flows. The background turbulence kinetic energy kQ is defined as a frac- 
tion of the mean flow kinetic energy; that is, in the present computation, 
kQ = (5.56 X I0“5)u|^ with u^y being the average velocity of the airstream. 
The turbulence kinetic energy due to the shear stress depends on the local 
maximum velocity gradient 3u/3s (3u/3s is the greater value of 3u/3y and 

3u/3z) and a mixing length Ijg. For flow near the duct wall, the mixing length 


19 



in the conventional turbulent boundary layer was used; near the jet, the mixing 
length in the axisymmetric mixing layer was assumed. Thus, 
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(16) 


In the present coordinate system (fig. 12), for flow near the duct wall, 
z ^ 0.01905 - 0. 096tj/0. 435 with 65 = 0.007 m being the wall boundary- layer 
thickness. In this region, Zjp is (ref. 11) 
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with 6 jjj = 0.005 m being the mixing-layer thickness. Similarly, the turbulence 
dissipation rate and the turbulent viscosity were assumed to be 
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No-flux boundary conditions were specified at the three symmetric bound- 
aries; on the duct wall, the velocity was zero, the temperature was constant, 
and heat transfer was allowed. 

In the finite-difference computation, nonuniform grid spacings were 
chosen; 11 x 30 grid points were used in the y-z plane. The computation was 
performed from the initial station (x = 0.0421 m downstream from the strut) 
to a downstream station (x = 0.7369 m) with the empirical constants - 1.44, 
C 2 = 1.92, Cf) = 0.09, Iciminar Prandtl number Npj-^£ = 0.7, and turbulent 
Prandtl numbers Npr,t,q =1*0 for q = u,v,w,k, Npr,t,e = 1*3» 

**Pr,t,H “ 0.9, and Np^^t,f =0.7. A total 690 variable-sized streamwise 
forward steps were taken. Total computer time was 261 sec on the Control Data 
CYBER 175 computer system; the total computer storage used was 774000 words. 

The accuracy of the computation was checked by the conservation of mass and 
mass flow. The mass source in each control volume was kept as small as 
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2 X 10~9 percent of the total mass flow, while the total mass-flow imbalance 
across the cross-stream plane at each station was kept under 2 x 10 ~^ percent 
of the total mass flow. Therefore, the total mass flow of the helium and air 
was virtually unchanged along the entire length of the duct during the 
computation. 


Finite-Element Program 

The finite-element program (COMOC) , has also been applied to the helium- 
air mixing experiment. Because of computer storage limitations and long run 
times when using COMOC, only one-quarter of the duct cross-sectional area was 
considered. Results for the remaining three quadrants associated with the jet 
were inferred by symmetry. The additional symmetry boundary condition was 
approximated by the initial data. There was little change in any of the depen- 
dent variables moving transversely between y = 0 and y = -0.00254 m. 

Initial data requirements for COMOC include the three velocity components, 
temperature, and species mass fraction profiles, and an initial estimate of 
streamwise pressures. The streamwise velocity component, temperature, and mass 
fractions were identical with those input into the finite-difference code. 

Small values of the v and w velocity components on the order of 0.1 percent 
of the streamwise velocity were assumed, consistent with continuity, to provide 
a gradient source on the y- and z-momentum equations as no transverse or lateral 
pressure gradient sources existed. (In the finite-difference program, such 
sources existed.) No effect was observed on the final downstream results from 
this assumption. The streamwise pressure gradient was initially assumed to be 
zero, but as disoussed previously, the pressure algorithm rapidly generated a 
gradient to match computational and actual cross-sectional areas. 

Turbulence closure was developed according to mixing-length theory. In 
this case, the turbulent viscosity was given by (ref. 2 ), 
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where C)^ = 0.4, X = 0.09, 65 is the boundary-layer thickness on the wall, 

and D is the Van Driest damping factor given by 
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vhere 


u(0. 01905 - z) 

z = (21b) 

V 



and T„ is the skin friction on the wallr the wall density, and v the 

kinematic viscosity. The effective Prandtl numbers for enthalpy and mass 
diffusion were, respectively, 1.0 and 0.4. 

Boundary conditions for the three symmetry boundaries were chosen to 
require that the normal gradient of each of the dependent variables across 
each of these boundaries vanish. Along the wall boundary, a gradient boundary 
condition was applied to the streamwise velocity component just off the wall 
at the outer edge of the sublayer boundary. This condition required that the 
normal derivative of u be equal at the sublayer boundary to the ratio of the 
local turbulent shear stress to local turbulent viscosity. Turbulent shear 
stress was found from a modified relation developed by Patankar and Spalding 
(ref. 13) through the application of wall functions within the sublayer. Thus, 
the "gradient boundary condition" simply inverted the basic definition of wall 
shear and solved for the velocity gradient given shear and viscosity. Although 
the program can carry the computation directly to the wall, a large number of 
elements were required to compute the large normal gradients encountered in the 
sublayer. Since the flow field is conveniently predicted in the sublayer by 
empirical means, such a computation to the wall was unjustified. The v and 
w velocity components were required to satisfy a no-slip condition at the wall; 
that is, they vanish "near" the wall. The normal gradients of both enthalpy 
and helium mass fraction were required to vanish at the wall, consistent with 
a no-flux condition there. 

The finite-element program also uses a nonuniform grid to discretize the 
cross-stream (y-z) plane. This plane is spanned by a 6 x 16 node mesh and 150 
triangular finite elements. To march the solution from the initial station 
(x = 0.0421 m) to the final downstream station (x = 0.7369 m) required 1326 sec 
of computer time on the Control Data CYBER 175 computer system with a storage 
of 2650000 words. 


COMPARISON OF NUMERICAL CALCULATIONS AND EXPERIMENTAL DATA 

Downstream computations (x = 0.7369 m) from both computer programs are 
compared with the experimental data in figure 14. Results for streamwise 
velocity, temperature, and helium mass fraction are plotted from the center 
plane (z = 0) to the duct wall (z = 0.01905 m) at five stations along the 
y-axis. Because the finite-element program calculated only one-quarter of the 
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center jet region, comparisons for this program are made only at three stations 
in the positive y-direction (figs. 14(a), 14(b), and 14(c)). 


Finite-Difference Program Results 

Comparisons of the finite-difference (SHIP) results with experimental data 
(fig. 14) generally show excellent agreement. The computed streamwise velocity 
and static temperature give excellent quantitative and qualitative agreement 
with the experimental data at all five stations along the y-axis. The maximum 
differences between computational and experimental data are about 6.3 percent 
for the velocity and 6 percent for the temperature. These maximum differences 
seem to occur at the measurement points near the duct wall. Several explana- 
tions for this discrepancy are possible. The experimentally measured values 
near the duct wall may be perturbed by the end effects of the duct which may 
propagate upstream through the subsonic wall boundary layer. Also, data reduc- 
tion assumed constant total temperature, and this assumption may introduce 
small errors. On the other hand, in the finite-difference computation, the 
deficiency near the wall may be due to both the approximate wall boundary con- 
dition and the coarse grid spacings selected near the wall. Despite these 
approximations and uncertainties, the comparisons near the wall are quite 
satisfactory. 

The computed helium mass concentrations are also in good agreement with 
the experimental data except at the two outermost stations (y = ±0.0127 m) . 

The disagreements at these stations may be attributed to several sources. 

First, the downstream measured helium mass flow is experimentally 11 percent 
higher than the upstream measured helium mass flow; however, the computed 
helium mass flow is conserved relative to the initial (upstream station) mass- 
flow data. Second, compared with the experimental data, the mixing of helium 
with air in the y-direction is relatively slow for the computation; this may 
be due to an improper estimate of the initial turbulence kinetic energy, dis- 
sipation rate, and (three-dimensional) mixing length. In addition, the flow 
field of the center jet mixing region was assumed to be independent of all 
other jets in the computation. Practically, there are interactions between two 
neighboring jet regions, particularly between the central jet and the cross-flow 
jet. These interactions may have caused nonsymmetric distribution of helium 
mass concentration at the stations in the positive and negative y-directions. 
Hence, the helium may be diffused in and out from the boundaries with neighbor- 
ing regions. Moreover, the absolute values of the helium mass fraction are 
actually quite small, on the order of less than 0.03. The accuracies of the 
measurement and the data reduction may limit such a quantitative comparison. 


Finite-Element Program Results 

Comparisons of the finite-element results with experimental data 
(figs. 14(a), 14(b), and 14(c)) also show generally good agreement. The static 
temperatures are in excellent agreement with the experiment, having a maximum 
difference of only 6.0 percent. The computed streamwise velocities also agree 
well with the experimental measurements, differing by a maximum of 8.6 percent. 
Values of the computed helium mass fraction are in good agreement with the 
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experimental data at the first two lateral stations (figs. 14(a) and 14(b))/ 
but agreement is poor at the outermost station (fig. 14(c)). In addition to 
the possible sources of disagreement discussed for the finite-difference 
results, the absence in the finite-element program of a transverse pressure 
gradient to convect helium away from the jet center line may be the major prob- 
lem. This deficiency would explain the need of a Prandtl number for mass dif- 
fusion (Schmidt number) which is lower than normally observed. Additionally, 
a coarse element discretization was necessary to maintain reasonable computer 
time requirements. As a result, the mass flow rates of helium and air were not 
conserved along the duct, and an overall computational loss in mass flow rate 
of 15 percent was detected. 


Comparison of Program Results 

In view of such a complicated three-dimensional flow field, the predictions 
from both computer programs are considered reasonably good in comparison with 
the experimental results. Comparisons between the two numerical computations, 
however, indicate that the finite-difference program provides better overall 
agreement with experimental data and requires less computer time than the 
finite-element program. 

The superior computational accuracy and run time of the finite-difference 
program are due primarily to the efficient implicit integrator used to solve 
the governing equations describing the flow field. This feature allows a very 
fine nodal discretization of the computational domain, permitting an accurate 
resolution of gradients within the flow field and requiring reasonable computer 
run times. The finite-element code currently uses a less efficient explicit 
predictor -corrector scheme. Recent improvement of the integration algorithm in 
the finite-element program has, however, resulted in a significant improvement 
in integration speed and run time. Better integration efficiency also allows 
the effective use of the two-equation turbulence model by the finite-difference 
program. The finite-element program is limited by its integrator to an alge- 
braic turbulence model. Finally, the finite-difference program calculates a 
three-dimensional pressure field, whereas the finite-element program calculates 
only an axial pressure field. Both transverse pressure gradients are important, 
however, in mixing flows. Therefore, a three-dimensional pressure algorithm is 
currently being added to the finite-element program. 


Comparison of Computational and Experimental Mixing Efficiencies 

In figure 15, the mixing efficiencies calculated along the duct by the 
two computer programs are presented. The mixing efficiencies at the upstream 
(x = 0.0421 m) and the downstream (x = 0.7369 m) measurement stations were also 
calculated directly from the experimental surveys and are plotted for compari- 
son. Mixing efficiency is a measure of the completeness of mixing between the 
injectant and the surrounding airstream amd is defined at any streamwise sta- 
tion as the fraction of fuel (injectant) that would react if complete reaction 
occurred without further mixing (ref. 15) . TO produce a consistent definition 
in the present calculation of mixing efficiencies, the helium mixed with the 
airstream was treated as if it were hydrogen, mixed and reacted with the same 
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airstream. The change in mixing efficiency along the duct is very useful for 
practical combustor design and helps to determine the combustor size required 
to achieve desired combustion per formetnce. In experiments, the mixing effi- 
ciency can only be obtained at relatively few stations where the flow fields 
are surveyed. In numerical computations, however, the mixing efficiency can be 
conveniently calculated at each station along the duct. Such calculations are 
extremely valuable to supplement the traditional approach based purely on 
exper imentation . 

In the comparison of the computational and experimental mixing efficien- 
cies (fig. 15), results again generally agree. However, the computational mix- 
ing efficiencies are higher than the experimental value at the initial station 
(x = 0.0421 m) . These differences are attributable to both physical and numer- 
ical reasons. The experimental mixing efficiency was integrated over the entire 
region around the central jet by using 819 integration points, while the finite- 
difference result was integrated over the upper half of the region by using 
330 nodes and the finite-element result used 96 nodes (150 finite elements) over 
one-quarter of the region. Because of strong gradients of helium mass fraction 
and velocity at the initial station, the use of a different number of integra- 
tion points naturally results in different calculated mixing efficiencies. 
Furthermore, the flow field in the central jet region is not exactly symmetric 
with respect to the symmetry boundaries assumed in the computations; hence, 
mixing efficiencies calculated in different regions will also be different. At 
the downstream station (x = 0.7369 m) , the mixing efficiencies of both computa- 
tions are in good agreement with the calculation made directly from the experi- 
mental survey data. 


CONCLUDING REMARKS 

Finite-difference and finite-element computer programs have been applied 
to a three-dimensional nonreacting turbulent mixing experiment. Relative accu- 
racies and efficiencies of the programs have been assessed by comparing their 
results with the experiment and also by comparing the results of the programs 
themselves. 

In general, both computer programs predict the experimental results and 
data trends reasonably well. Their predictions were based on the same set of 
initial experimental data for streamwise velocity, temperature, pressure, and 
concentration, but the programs used different turbulence models, governing 
equations, and pressure algorithms. Compared with experimental data at a down- 
stream station, the predictions by the finite-difference program are generally 
more accurate than those of the finite-element program. Additionally, the 
finite-difference program is more efficient in both computer storage and com- 
puting time than the finite-element program. The superior accuracy of the 
finite-difference program is due primarily to its efficient integrator, which 
allows a very fine discretization of the flow field and the application of a 
more appropriate turbulence closure model. Additionally, the finite-difference 
program calculates transverse pressure gradients in the duct, whereas the 
finite-element program is limited to a one-dimensional streamwise pressure 
field. Both limitations degrade the accuracy of the finite-element program. 
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and work Is currently underway to upgrade both the integrator amd the pressure 
algorithm. 

At present, the finite-difference program is preferred for numerical com- 
bustor cinalysis and design. Development of the finite-element program has, 
however, been underway for significantly less time. Developed to their full 
potential, both computer programs should provide powerful analysis tools for 
scramjet engine design. 


Langley Research Center 

National Aeronautics eind Space Administration 
Hampton, VA 23665 
February 16, 1978 
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APPENDIX 


ASSEMBLED SOLUTION ALGORITHM 

In this appendix, the assembled finite-element solution algorithm is 
developed. To simplify this development, the following notation, which corre- 
sponds more closely to that of reference 2, is used: 

nj^ unit vector normal to kth coordinate direction 

streamwise coordinate x 

*2' *3 cross-section coordinate corresponding to y and z, respectively 
Subscripts: 

i coordinate index ranging from 1 to 3 

k coordinate index taking on values of 2 and 3 

The summation convention is also used. 

In this notation, equation (11) becomes 

\ {<MXk)}L(q*) dV- X(p {(Mxk)}JMqm) dS = {0} (Al) 

'^Rm ^aR„n3k 

and equations (7) become 


L(q) = ic 


8xi 


K(q) 


3xl 


/ 3q \ 

+ flq,— ,Xij - g(q,xi) = 0 






(A2) 


S,(q) = a(^) q(x 2 ,Xk) + a^^) K(q) q(X 2 ,Xk) nk - a^^^ = 0 

9*k ^ 


Now, substituting equations (A2) for L and S, into equation (Al) gives 


\ {<t>(Xk)}K 


9xk 


K(qii) 


8xk 


dV + \ {<t> (xk) } 


'*ij “ 9(qm'*l) 


dV 


- X® {(|>(Xk)} 

^3R„,n3R 


( 1 ) 


,( 2 ) 


- 


% 9m(*l'*i) ®ra - 9m(*l»*i)*'k ~ ®m 

3xk 


dS = {0} 
(A3) 
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Ekjuation (A3) presents a potential problem as a consequence of the method of 
weighted residuals. The differential operators L and H contain deriva- 
tives as high as second order with the possibility of introducing discontinu- 
ities into the integrands at points in the flow field. The certainty of finite 
derivatives can be improved if the order of L and 5. can be reduced. 

The Green-Gauss theorem, given by 



(A4) 


provides a means for reducing the order of the operators. The first term of 
equation (A3) can be modified to fit the form of equation (A4) as follows: 


I 


{({)}< — 

Rm 




dV 


4 - 




3q, 


im 


9xr 


dV 


- K ) K dV 

3xk 




^ , 3 3<3m 3 / \9qm 

{<t)}K + ({(t>}K) 

3x|^ 3 x) 5 3xj^ \ / 3xjj 


H£ 


3{(j)} 9gm 

K 

Rm 9xk 


dV 


(A5) 


Note that the functional notation has been dropped to simplify the equations. 
By identifying $ with ({(|)}k) and with qj„, the Green-Gauss theorem can 
be applied to the first term of equation (A3) to give 


I. 


{4)}ic — 
Rm 3^k| 


3xk 


dV 


9% 


= icd^ {(|)}K nu dS - K r 

3X]j 


3{(|)} 9qm 


X 9xk 3xk 


dV 


(A6) 
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Substituting equation (A6) into equation (A3) then gives 


K J {<t>}K — nk dS 


9Rjn 


3xk 


-< r 




dV + y {<t.}[f* - g*] dV 


x(5 4 } 

•^3RnP3R 


( 1 ) 




qm + ai 


( 2 ), 


3q, 


m 


m 


3xk 


n 

•'k ~ 


dS = {0} 


(A7) 


where = f ^qm» 3qmy3*i'*i) 9m “ Note that the derivative 

order of equation (A7) has been lowered by 1 with respect to equation (A3) . 

(2) 

Now, since X is an arbitrary algebraic multiplier, choose AajJl equal to 
K for 3Rm n 9R nonvanishing (ref. 7) ; this choice causes the first term of 
equation (A7) to cancel with the second term of the fourth integral only where 
3Rjn and 3R correspond. This leaves 




3qm 

<4^ {({)}K 

9Rm 

3R„in3R=0 


dS - K J ^ K ^ dV + J W[f* - g^] 

^Rfl, 9*k 9*k *^R„, ^ -* 


dV 


- ... '♦>[ 


3Rnfi3R 


r.-,. (1)' * ^(3) 

3m qm “ 3m 


■] 


dS = {0} 


(A8) 


where 




and 




Equation (A8) provides the solution algorithm for each finite-element sub- 
domain. To establish the solution for the entire duct cross-sectional plane 
(the global solution), equation (A8) is first evaluated for each of the M 
finite elements (m = 1, 2, . . ., M) . For each element, three first-order 
ordinary differential equations result for each dependent variable. These 
three equations are now summed (assembled) over the M finite elements accord- 
ing to Boolean algebra (symbol u) to give a 3 x m equation system that will 
produce the global solution. Thus, the global solution is 
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U 




3qm 

{(j>}K njj dS 

3% 

31^3 R=0 


w(4-gm)<3v 

*^Rm 3*k ‘^Rj„ ' ' 


ic $ 4}(am^^ 

'^3RnT>3R 


qm ” I 


•) 


= { 0 } 


(A9) 


Now, in accordance with reference 14, since {(j)} is a linear function, the 
surface contributions resulting from the first integral of equation (A9) can 
be made to cancel in pairs upon assembly. This result is physically meaningful 
since the evaluation of the surface integral with outward pointing normal about 
a common side of adjacent element pairs would give equal and opposite values 
that would cancel when their algebraic sum was taken. Then, equation (A9) 
becomes 


u 


■"I 


3{<t)} 3q: 


Pro 3xk 


{<t)}(fm - 9m) dV 
®m 


<C5 


ta qm “ “jn 


') 


dS 


= {0} 


(AlO) 


After the notation is changed to that of the main text, equation (AlO) is 
identical to equation (12) and is the basis of the solution algorithm for all 
dependent variables in COMOC. 
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TABLE II.- EXPERIMENTAL DATA AT UPSTREAM STATION (x = 0.0421 m) 


y» 


u, 

P, 

T, 

fHe» 

m 

m 

m/sec 

MPa 

K 

percent 

-0.0127 

0.0173 

436 

0.1172 

201 

0 


.0127 

497 

.1145 

171 

.0484 


.0064 

427 

.1124 

203 

.0971 


0 

460 

.1131 

189 

0 


-.0058 

457 

.1124 

191 

0 


-.0119 

451 

.1145 

193 

.0318 


-.0160 

280 

.1172 

256 

0 

-0.0102 

0.0173 

424 

0.1172 

206 

0 


.0127 

500 

.1145 

171 

0 


.0064 

433 

.1124 

202 

.0762 


0 

485 

.1131 

182 

.873 


-.0058 

460 

.1124 

189 

.00138 


-.0119 

451 

.1145 

192 

0 


-.0160 

295 

.1172 

251 

0 

-0.0076 

0.0173 

427 

0.1172 

204 

0 


.0127 

494 

.1151 

174 

0 


.0064 

436 

.1124 

201 

.0804 


0 

512 

.1124 

186 

4.67 


-.0058 

457 

.1124 

191 

.0443 


-.0119 

451 

.1151 

193 

0 


-.0160 

311 

.1172 

246 

0 

-0.0051 

0.0173 

457 

0.1172 

191 

0 


.0127 

494 

.1151 

173 

0 


.0064 

436 

.1124 

203 

.972 


0 

750 

.1124 

163 

26.7 


-.0058 

454 

.1124 

192 

.0166 


-.0119 

451 

.1151 

193 

.00138 


-.0160 

305 

.1172 

248 

0 

-0.0025 

0.0173 

469 

0.1172 

185 

0 


.0127 

494 

.1151 

174 

.289 


.0064 

472 

.1124 

204 

5.27 


0 

1152 

.1117 

143 

80.2 


- . 0058 

451 

.1124 

195 

.337 


-.0119 

445 

.1151 

196 

0 


-.0160 

348 

.1172 

234 

0 
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TABLE II.- Concluded 


y. 

m 

z. 

m 

u, 

m/sec 

P. 

MPa 

T, 

K 

f'He* 

percent 

0 

0.0173 

469 


184 

0 


.0127 

491 

■ms 

175 

0 


.006H 

518 

.1124 

199 

9.69 


0 

1140 

.1103 

143 

77.7 


-.0058 

451 

.1124 

195 

.367 


-.0119 

445 

.1145 

196 

0 


-.0160 

314 

.1165 

246 

0 

0.0025 

0.0173 

463 

0.1151 

188 

0 


.0127 

500 

.1131 

170 

0 


.0064 

445 

.1117 

199 

.760 


0 

686 

.1089 

170 

20.9 


-.0058 

457 

.1117 

190 

.0429 


-.0119 

454 

.1131 

192 

0 


-.0160 

320 

.1151 

244 

0 

0,0051 

0.0173 

454 

0.1103 

192 

0 


.0127 

503 

.1117 

169 

0 


.0064 

439 

.1110 

199 

.0762 


0 

500 

.1069 

179 

1 .80 


-.0058 

460 

.1110 

189 

.00276 


-.0119 

457 

.1117 

191 

0 


-.0160 

323 

.1103 

243 

0 

0.0076 

0.1173 

433 

0.1089 

201 

0 


.0127 

506 

.1117 

168 

0 


.0064 

439 

.1103 

199 

.0734 


0 

482 

.1055 

179 

.214 


-.0058 

466 

.1103 

186 

.00276 


-.0119 

460 

.1117 

189 

0 


-.0160 

344 

.1089 

236 

0 


0 

0 
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TABLE III.- EXPERIMENTAL DATA AT DOWNSTREAM STATION (x = 0.7369 m) 




0.0173 

.0127 

.0064 

0 

-.0058 

-.0119 

-.0160 

0.0173 

.0127 

.0064 

0 

-.0058 

-.0119 

-.0160 


0.0173 

.0127 

.0064 

0 

-.0058 

-.0119 

-.0160 

0.0173 

.0127 

.0064 

0 

-.0058 

-.0119 

-.0160 


0.0173 

.0127 

.0064 

0 

-.0058 

-.0119 

-.0160 


m/sec 


0.1793 

.1703 

.1648 

.1669 

.1600 

.1613 

.1641 











TABLE IV.- COMPARISON OF THE TWO COMPUTATIONS BASED ON SHIP AND COMOC 



SHIP 

COMOC 

Governing equations 

Parabolic type of Navier- 
Stokes equations 

Parabolic type of Navier- 
Stokes equations with 
additional boundary- 
layer-type simplifications 

Pressure gradient 

Streamwise and transverse 

Streamwise only 

Turbulence model 

"k-e" two-equation model 

Mixing-length model 

Reaction model 

Equilibrium reactions 

Complete reaction 

Numerical scheme 

Finite difference 

Finite element 

Region of computation 
(fig. 12) 

ABCDEFA 

ABCFA 

Number of grid points 

11 X 30 

6 X 16 

Input initial flow 
field data 

u, P, T, fHe 
(v = w = 0) 

u, p, T, f°^ (v, w estimated 
from continuity equation) 

Number of streamwise 
stations computed 

690 steps 

351 steps 

Total computer 
storage 

774000 

2650000 

Computing time on 
Control Data 
CYBER 175 

261 sec 

1326 sec 
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All dimensions are in meters. 



Projection of injectors 


Upper wall 



Figure 6.- Helium mass fraction contours of duct cross section at x r 0.7369 m. 









Lower wall 


: = 0.7369 ni. nifjg o 0 
Figure 8.- Concluded. 
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Lower wall 


(a) X = 0.0421 m. 

Figure 9.- Streamwise velocity contours of center jet mixing region. 

Velocity is given in m/sec. 








Lower wall 


(a) X B 0.0^121 m. 

Figure 10.- Static pressure contours of center jet mixing region 
Static pressure is given in mPa. 






Upper wall 



Lower wall 

(b) X = 0.7369 m. 

Figure 10 .- Concluded 
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Lower wall 


(a) X = 0.0M21 m. 

Figure 11.- Static temperature contours of center jet mixing region 
Static temperature is given in K. 
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(b) X B 0.7369 m. 


Figure 11 .- Concluded 
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(a) y = -0.0076 m. 

Figure 13.- Experimental flow variables at initial station (x = 0.0421 m) . 

Solid lines indicate interpolation. 
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Figure 13.- Continued. 
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(c) y = -0.0025 m. 


Figure 13.- Continued. 
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Figure 13.- Continued. 
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Figure 13.- Continued 
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(f) y = 0.0051 m. 
Figure 13.- Continued. 









O □ Experimental data 
Finite difference 



-2 -2 

f^e> percent, T x lo , K, u x 10 , m/sec 


(b) y = 0.0063 m. 
Figure 14.- Continued. 



O O t> Experimental data 
Finite difference 



f^g, percent, T x lo'^, K, u x io“^, m/sec 


(c) y H 0.0127 m. 
Figure 14.- Continued. 







O □ P> Experimental data 
Finite difference 



percent, T x io~^, K, u x lo’^, m/sec 

(e) y = -0.0127 m. 

Figure 14.- Concluded. 
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